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ABSTRACT 

Direct imaging of exoplanets is limited by bright quasi-static speckles in the point spread function (PSF) of 
the central star. This limitation can be reduced by subtraction of reference PSF images. We have developed 
an algorithm to construct an optimized reference PSF image from a set of reference images. This image is 
built as a linear combination of the reference images available and the coefficients of the combination are 
optimized inside multiple subsections of the image independently to minimize the residual noise within each 
subsection. The algorithm developed can be used with many high-contrast imaging observing strategies relying 
on PSF subtraction, such as angular differential imaging (ADI), roll subtraction, spectral differential imaging, 
reference star observations, etc. The performance of the algorithm is demonstrated for ADI data. It is shown 
that for this type of data the new algorithm provides a gain in sensitivity by up to a factor 3 at small separation 
over the algorithm used in Marois et al. (2006). 

Subject headings: Instrumentation: adaptive optics — planetary systems — stars: imaging — techniques: 
image processing — techniques: high angular resolution 
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1. INTRODUCTION 

Direct imaging of exoplanets, circumstellar disks, jets, 
winds or other structures around stars is difficult due to the an- 
gular proximity of the star and the very large luminosity ratios 
involved. Current attempts, both from the ground with adap- 
tive optics (AO) and from space, are limited by a swarm of 
bright quasi-static speckles that completely mask out the faint 
planets or structures that are sought after (Schneider & Silver- 
stone 2003; Biller et al. 2004; Marois et al. 2005; Masciadri 
et al. 2005). These speckles are caused mainly by imperfec- 
tions in the optics and are long-lived, hence the "quasi-static" 
appellation. As the exposure time is increased, the quasi- 
static speckles add coherently and their intensity eventually 
becomes dominant over signals that add incoherently, such as 
sky or read noise and general (non-static) speckles. 

This problem is more important closer to the star, as the rel- 
ative speckle intensity is higher there, and the size of the re- 
gion in which the noise is dominated by quasi-static speckles 
will depend upon the exposure time, the sky and read noise 
levels, the telescope and camera used, the target brightness, 
etc. For example, the observations of Masciadri et al. (2005) 
obtained at the Very Large Telescope are limited by speckles 
only at subarcsecond separations, while in the search for plan- 
ets on wide orbits that we are currently carrying at the Gemini 
telescope (D. Lafreniere et al., in preparation), the observa- 
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tions, which use longer individual exposure times (30 s), are 
typically limited by quasi-static speckle noise out to separa- 
tions of ~5"-10". Quasi-static speckles even dominate the 
noise at separations well past 10" in the observations of the 
bright star Vega obtained by Marois et al. (2006) at the Gem- 
ini telescope; this also appears to be the case for similar ob- 
servations obtained on the Keck and Palomar 5-m telescopes 
(Macintosh et al. 2003; Metchev et al. 2003). At a given an- 
gular separation, no gain in contrast is achieved by increasing 
the exposure time once the noise is dominated by quasi-static 
speckles. 

When this regime is reached, it is possible to subtract the 
quasi-static speckles by using reference point spread function 
(PSF) images. A reference PSF image is any image whose 
subtraction from the target image would reduce the signal 
from the speckles while preserving that of the object sought 
after. For example, reference PSF images can be obtained 
from observations of reference stars, or from observations of 
the target itself obtained at different field of view orienta- 
tions (e.g. Schneider & Silverstone 2003), wavelengths (e.g. 
Racine et al. 1999), or polarizations (e.g. Kuhn et al. 2001). 

Obtaining a reference PSF image highly correlated with the 
target image is a difficult task because even though quasi- 
static speckles are long lived, they still vary with time due to 
temperature or pressure changes, mechanical flexures, guid- 
ing errors or other phenomena (Marois et al. 2005, 2006). On 
the other hand, even when a reference PSF image is acquired 
simultaneously with the target image at other wavelengths 
or polarizations, differential aberrations within the camera 
decorrelate the PSFs (Marois et al. 2005; Lenzen et al. 2004). 
Thus, when trying to subtract speckles one must always work 
with slightly decorrelated reference PSF images and the spe- 
cific way in which the available data are used to perform the 
subtraction may have a significant impact on the speckle noise 
attenuation achieved. This paper presents a way of combin- 
ing reference PSF images to optimize the noise attenuation. 
In particular the algorithm is applied to angular differential 
imaging (ADI) (Marois et al. 2006), which is currently one of 
the most efficient quasi-static speckle suppression technique 
for ground-based observations. Although emphasis is given 
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to point source detection throughout the paper, the algorithm 
can be optimized to search for any other structure in the close 
vicinity of a star. 

The new reference PSF construction algorithm is presented 
in §2. Then, a review of ADI and the algorithm used by 
Marois et al. (2006) is presented in §3. In §4, the new al- 
gorithm is applied to ADI and its performance is presented. 
The possibility of using this algorithm with other observing 
strategies is finally discussed in §5. 

2. REFERENCE PSF CONSTRUCTION BY LOCALLY 
OPTIMIZED COMBINATION OF IMAGES 

Consider a single target image, from which speckles are 
to be subtracted, and suppose that N reference PSF images 
are available for this purpose. The heart of the algorithm de- 
scribed here is to divide the target image into subsections and 
to obtain, for each subsection independently, a linear combi- 
nation of the reference images whose subtraction from the tar- 
get image will minimize the noise. By optimizing the weights 
given to the N available reference PSF images according to 
the residual noise obtained, this approach produces a repre- 
sentation of the target PSF image that is better than any pre- 
defined combination of the reference PSF images. Further, it 
is advantageous to optimize the coefficients of the linear com- 
bination for subsections of the image because the correlation 
between the target and the reference PSF images generally 
varies with position within the target image. We refer to the 
algorithm described here as "locally optimized combination 
of images", or LOCI. 

The coefficients used for subtraction of the speckles within 
subsection S T of the target image are determined by a mini- 
mization of the noise within a generally larger, so-called op- 
timization, subsection T , which encompasses S T . The cor- 
responding optimization subsections in the reference PSF im- 
ages are denoted O n , n= 1,. . . ,N. 

Ideally, to achieve the optimal noise attenuation everywhere 
in the target image, one would want to optimize the coeffi- 
cients for subsections S T that are as small as possible, ulti- 
mately consisting of a single pixel. In practice, to avoid a 
computationally prohibitive repetition of the algorithm, one 
uses subsections that contain many pixels, within which the 
same linear combination of reference images is used. 

While the size of the subsection S T is limited by computa- 
tion resources, the size of T is determined by the need to pre- 
serve the signal from any point source sought after. From the 
point of view of the algorithm described below, a point source 
in T is a residual that it tries to minimize and will partially 
subtract. The amount of partial subtraction depends upon the 
fractional area of T that is occupied by the point source. So, 
even though smaller optimization subsections lead to a better 
noise attenuation, they also lead to a larger subtraction of the 
signal of the point sources sought after. Thus the size of T 
must be properly determined and the amount of partial sub- 
traction of point sources must be well characterized. The area 
A of the optimization subsection is determined by the param- 
eter Na through the expression 

A=N A n(j) (1) 

where W is the full-width-at-half-maximum (FWHM) of the 
PSF; N A thus corresponds to the number of "PSF cores" that 
fit in the optimization subsection. 
If the set of reference PSF images contains target images, it 



is necessary to construct the optimized PSF to be subtracted 
from a given subsection S T by using only the subset of these 
images in which a companion point source appearing in S T 
would be displaced by at least a distance <5 m j n or would have 
an intensity smaller by at least a factor a with respect to its 
position or intensity in the image from which speckles are to 
be subtracted. In other words, this subset includes all refer- 
ence PSF images of index k£K, where 

K = {k€[l,N] : |r fc -r r | >5 mfa V f k /f T < a}, (2) 

being any field position in the subtraction subsection of 
the target image and r* the corresponding position in image 
k, while fk/fr is the intensity ratio in those images of any 
companion sought after. If the set of reference PSF images 
does not contain target images, then K = {1, . . . ,N}. The pa- 
rameters S m i n and a, when applicable, affect both the speckle 
noise attenuation and the amount of partial subtraction of 
point sources, similarly to A^. The best values to use, which 
depend on the type of data being analyzed and the level of 
correlation between the target and reference images, may be 
determined from a comparison of the results obtained with 
different values, see §4. 1 . For the remainder of the section, it 
is assumed that values for N A , S m in, and a have been selected 
by the user. 

The reference PSF for the optimization subsection is then 
constructed according to 

(f = Y J C k O k , (3) 

keK 

where the coefficients c k are to be determined by the algo- 
rithm. They are computed by minimizing the sum of the 
squared residuals of the subtraction of R from T , which 
is given by 

ai = ^ mi (Ol-Of) 2 = ^ mi (oJ-^c k O?j , (4) 

where i denotes a pixel in the optimization subsection and m 
is a binary mask that may be used to ignore some pixels. The 
quantity to minimize is a sum and can be biased by cosmic ray 
hits or bad pixels if they have not been properly corrected or 
filtered before the algorithm is used. When bad pixels remain 
in the image, the bias can be completely remedied by setting 
the mask m to zero for these pixels. Generally, the fraction of 
pixels affected is small and their exclusion from the computa- 
tion of the residuals has practically no impact on the solution 
found. The minimum of a 2 occurs when all its partial deriva- 
tives with respect to the coefficients c k are equal to zero, i.e. 
when 

U=E- 2 ^/(°r-E^)= ' v ^*- w 

Reversing the summation order and rearranging the terms we 
find 



^\^m i OlO k yY j m i O j i O T i , VjeK. (6) 

This is a simple system of linear equations of the form Ax = b 
where 
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Afl^miOjot, x k = c\ bj^miOjoJ. (7) 

/ i 

Solving this system gives the coefficients c k needed to con- 
struct the optimized reference PSF image for the subsection 
S T . By construction, assuming that all the O k are linearly in- 
dependent, the matrix A is always invertible. Thus, the system 
always has a unique solution, meaning that for the given op- 
timization subsection and set K the solution found leads to 
an absolute minimum of the residuals. Finally, using the set 
of optimized coefficients, the optimized reference PSF image 
subsection to be subtracted from S T is constructed as 

S R = J2c k S k , (8) 

where S k denotes the corresponding subtraction subsection in 
the reference PSF image k. 

3. REVIEW OF ANGULAR DIFFERENTIAL IMAGING 

The ADI technique, as detailed in Marois et al. (2006), con- 
sists in acquiring a sequence of many exposures of the target 
using an altitude/azimuth telescope with the instrument ro- 
tator turned off (at the Cassegrain focus) or adjusted (at the 
Nasmyth focus) to keep the instrument and telescope optics 
aligned. This is a very stable configuration and ensures a high 
correlation of the sequence of PSF images. This setup also 
causes a rotation of the field of view (FOV) during the se- 
quence. For each target image in such a sequence, it is pos- 
sible to build a reference image from other target images in 
which any companion would be sufficiently displaced due to 
FOV rotation. After subtraction of the reference image, the 
residual images are rotated to align their FOV and co-added. 
Because of the rotation, the PSF residual speckle noise is aver- 
aged incoherently, ensuring an ever improving detection limit 
with increasing exposure time. 

In building a reference image, a compromise has to be 
reached between the quasi-static speckle noise correlation, 
which is highest for the shortest time delays between images, 
as shown in Figure 2 of Marois et al. (2006), and the need to 
ensure a sufficient companion displacement. The minimum 
time delay T m j n between an image and the ones which can be 
used as references decreases as the inverse of the angular sep- 
aration. Accordingly, it is possible to use images more closely 
separated in time to build the reference image at larger angular 
separations. 

In the speckle subtraction algorithm used in Marois et al. 
(2006) (see their §5.2 and their Table 2), the first step is to 
subtract the median of all the images from each individual 
image. Each target image is then broken into many annuli to 
reflect the dependence of r m i n on the distance from the center 
of the PSF. A reference image is obtained within each annulus 
by median combining the four images obtained closest in time 
but at least T m [ n from the target image. The intensity of this 
reference image is also scaled to minimize the noise after ref- 
erence subtraction. All the resulting images are then rotated 
to align their FOV and a median is taken over them. 

4. APPLICATION OF THE LOCI ALGORITHM TO 
ANGULAR DIFFERENTIAL IMAGING DATA 

4. 1 . Definition of the arbitrary parameters specific to ADI 

As mentioned in §2, some parameters must be chosen by 
the user before the algorithm is used. For ADI data, as images 




FIG. 1. — Example of subtraction (shaded in grey) and optimization (de- 
limited by thick lines) subsections for ADI using the procedure of §4.1. The 
left and right panels show the subtraction and optimization subsections for 
the 1 st and 13 th subtraction annuli respectively. In the right panel, the first 
12 subtraction annuli (of width dr) are marked by thin lines; dr increases 
with separation in this specific example. The central circle (cross-hatched) 
represents the saturated region. 

are generally acquired in a single bandpass, the parameter a 
does not apply. On the other hand, the area and shape of the 
subtraction and optimization subsections must be defined as 
well as the minimum displacement <5 m j n between sources in 
the target and reference images. 

The dependence of T m j n on angular separation suggests 
the use of annular geometry for the subtraction subsections. 
These subsections are obtained by further dividing the annuli 
azimuthally to reduce their spatial extent, which enables a bet- 
ter fit of local PSF variations as explained above. Since r m i n is 
proportional to l/r, the set of images that can be used to con- 
struct a reference PSF changes rapidly with radius at small 
separation and it is best to use narrow subsections at small 
radii to ensure that the largest possible set of reference PSF 
images is used at any separation. The subtraction subsections 
S T are described by their inner radius r, mean angular position 
4>, radial width dr, and angular width A(f>. 

The optimization subsections T share the same inner ra- 
dius, mean angular position, and angular width as their corre- 
sponding subtraction subsection. As explained in §2, the area 
A of the optimization subsections has to be chosen to max- 
imize noise attenuation while minimizing point source sub- 
traction. For an annular subsection T of inner radius r, radial 
width Ar, and azimuthal width (r+ Ar/2)A0, 

A = Ar(r+Ar/2)Acj). (9) 

We define the ratio of the radial and azimuthal widths of the 
optimization subsections as 

^At/W ^ (10) 

Then, by Eq.(l), Ar and A(f> are uniquely determined by the 
parameters g and N A and the PSF width W: 

Ar=^TTgN A W 2 /4, (11) 

and 

The optimization subsections were chosen not to be centered 
radially on the subtraction subsections but to extend to larger 
radii because, in the optimization, the radial dependence of 
the PSF noise gives more weight to the inner pixels, i.e. to the 
pixels in S T . Figure 1 shows an example of subsections that 
can be used with this procedure. 
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For ADI data, the set of target images is the same as the 
set of reference PSF images. The subset of images that can 
be used as references for a given subsection S T of a given tar- 
get image depends upon the parameter S m i„ introduced in §2, 
whose value is set by the parameter^ through the expression 

6 mm = N s W + rd6 n , (13) 

where d6 n is the angle of FOV rotation that occurred during 
exposure n. The last term of the expression above represents 
the azimuthal smearing of an off-axis point source that oc- 
curs during an exposure due to FOV rotation. The parameter 
Ns represents the minimum gap allowed, in units of the PSF 
FWHM, between a source position in image n and the corre- 
sponding positions in the images used as references. 

The values of Na, g, dr and Ns that maximize the sensitivity 
to faint point sources will be determined in the next section 
using real data. 

4.2. Parameter determination 

Observations of the star HD97334b (GOV, H = 5) were 
used to determine the values of the algorithm parameters. 
These observations are part of the Gemini Deep Planet Sur- 
vey (GDPS, D. Lafreniere et al., in preparation), which is 
an ongoing direct imaging search for Jupiter mass planets on 
large orbits (> 40 AU) around young nearby stars (—100 Myr, 
<35 pc). This particular dataset consists in a sequence of 90 
30-s images in the CJ-U-short (1.58 /im, 6.5%) filter obtained 
with ALTAIR/NIRI at the Gemini North telescope (program 
GN-2005A-Q-16). The f/32 focal ratio of the camera and 
8-m primary mirror diameter lead to 0."022 pixel -1 . The im- 
ages are saturated inside a radius of — 0"7 from the PSF cen- 
ter. Short unsaturated exposures were acquired before and 
after the saturated sequence to calibrate photometry and de- 
tection limits. The corrected PSF FWHM was measured to 
be 3.4 pixels, or 0"074, and the Strehl ratio was -16%. The 
Cassegrain rotator was fixed during all observations. Basic 
image reduction and registering was done as in Marois et al. 
(2006). 

The same procedure was used for optimizing each of Ns, 
Na, g and dr. First, the unsaturated PSF image was used to 
introduce artificial point sources to the images at angular sep- 
arations in the range 50-300 pixels (27-160 X/D) in steps of 5 
pixels (2.75 X/D). Each artificial source was smeared accord- 
ing to its displacement during an integration, and its intensity 
was set so that its S/N would be —10 in the final residual com- 
bination. Next, a symmetric radial profile was subtracted from 
each image to remove the seeing halo. Then the subtraction 
algorithm was run on the sequence of images with a range 
of values for the parameter under consideration. Finally, the 
noise and the flux of each artificial point source in an aper- 
ture diameter of one FWHM were measured in the residual 
image. This process was repeated 50 times by placing the ar- 
tificial sources at different angular positions each time. The 
trial values for the optimization of each parameter are listed 
in Table 1 . For dr either a fixed value is used throughout the 
image or one that varies from 1.5 to 15 in units of the PSF 
FWHM. The optimal value of a parameter was determined re- 
cursively, with the values of the other parameters set first to 
the medians of the values listed in Table 1 and then to their 
most recently determined optimal value except for dr set at a 
fixed value of 1.5. The results are shown in Figures 2-5. 

As can be seen in Figure 2, the minimum spacing has little 
impact on the recovered flux at separations > 100 X/D, where 



TABLE 1 

Parameter values used for optimization 



Parameter 


Trial values 


Adopted value 


N S 


0.25,0.5,0.75, 1.0, 1.5,2.0 


0.5 




50, 100, 150, 300, 500 


300 


g 


0.5, 1.0, 2.0 


1.0 


dr 


1.5,3,6,9, 15, (1.5-15) a 


(1.5-15) a 



a dr equal to 1.5 for separations less than 60 \/D and 15 for 
larger separations. 
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FIG. 2. — Average residual intensity of the artificial point sources normal- 
ized to their initial intensity (top), and their S/N ratio (bottom) as a function 
of angular separation, for different values of Ng . The solid, dotted, dashed, 
dot-dashed, triple-dot-dashed, and long-dashed curves are respectively for 
N s = 0.25, 0.5, 0.75, 1.0, 1.5, and 2.0. 
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FIG. 3. — Average normalized residual intensity (top) and S/N ratio (bot- 
tom) as a function of angular separation for different values of Na ■ The solid, 
dotted, dashed, dot-dashed, and triple-dot-dashed curves are respectively for 
N A = 50, 100, 150, 300, and 500. 

—80-90% of the flux is recovered independently of Ns- How- 
ever, at small separations the effect is important and signifi- 
cant loss in signal occurs, particularly for the smallest mini- 
mum displacements. This is because the fraction of images in 
the subset K for which the point source partially overlaps that 
in the target image is greater for smaller separations, where 
linear motion of the point source is slower. The best overall 
S/N is obtained with Ns = 0.5, for which the loss in the re- 
covered flux is more than compensated by the improvement 
in quasi-static speckle noise attenuation. 

Figure 3 shows that the residual signal of point sources is 
strongly dependent upon the size of the optimization subsec- 
tions, as expected from the discussion in §2. When Na is too 
small, the gain in attenuation is not sufficient to compensate 



PSF subtraction algorithm for high-contrast imaging 



5 




at , , , , , , 

40 60 80 100 120 140 160 

Angular separation {X/D) 



FIG. 4. — Average normalized residual intensity (top) and S/N ratio (bot- 
tom) as a function of angular separation for different values of g. The solid, 
dotted, and dashed curves are respectively for g = 0.5, 1, and 2. 



Angular separation (arcsec) 



6 



















! / /' ! 
\ i i \ 



40 60 80 100 120 140 160 

Angular separation {X/D) 



FIG. 5. — Average normalized residual intensity (top) and S/N ratio (bot- 
tom) as a function of angular separation for different values of dr. The solid, 
dotted, dashed, dot-dashed, triple-dot-dashed, and long-dashed curves are re- 
spectively for dr = 1.5, 3, 6, 9, 15, and dr varying with radius (see text). 

for the larger point source subtraction and lower S/N ratios are 
obtained, especially at large separations. On the other hand, 
when Na is too large, the quasi-static speckles are not sub- 
tracted as efficiently at small separations and lower S/N ratios 
result. A value of Na = 300 provides the best overall S/N ratio. 

The parameter g has little effect on the performance, see 
Figure 4, although for angular separations < 50 X/D regions 
more extended radially (g = 2) fare slightly better than regions 
more extended azimuthally. Nevertheless, we adopt g = 1 as 
the optimal value. 

Finally, Figure 5 shows that at small separations (< 
60 X/D), a dr > 6 leads to a lower S/N ratio because it poorly 
matches the evolution of T m j n with separation, as expected. 
Since a larger dr leads to a faster execution of the algorithm, 
because fewer subtraction subsections are required to cover 
the entire image, we use as the optimal value a dr equal to 1.5 
for separations less than 60 X/D and 15 for larger separations. 

The optimal parameter values may vary slightly from those 
found above for other sets of data depending on the telescope, 
instrument, seeing, FOV rotation rate, target brightness, etc. 
They are optimized here for a specific set of data only to il- 
lustrate the potential of the LOCI algorithm for ADI. For all 
computations that follow, the optimal values listed in Table 1 
are used. 

4.3. Point source photometry 

Since the algorithm reduces the flux of point sources sig- 
nificantly, especially at small separations, it is important to 
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FIG. 6. — Average normalized residual intensity (top) and ratio of the mea- 
sured dispersion of the residual intensity of sources over the residual noise 
(bottom). The solid, dotted, dashed, and dot-dashed lines are respectively for 
point source intensities yielding S/N of 3, 6, 10, and 25 in the final residual 
image. 

verify that the true flux can be recovered accurately and that 
the uncertainty on this value can be well determined. The 
algorithm was run on the sequence of images, with artificial 
companions of various intensities added at all angular separa- 
tions in the range 50-300 pixels (27-160 X/D) by steps of 5 
pixels (2.75 X/D). Four intensities were used, yielding S/N of 
3, 6, 10 and 25 in the final residual image. This process was 
repeated 50 times with the sources at different azimuthal po- 
sitions. The mean normalized residual source intensities and 
residual intensity dispersions over the 50 azimuthal positions 
were then computed and the results are shown in Figure 6. 

The top panel of this figure shows that the normalized resid- 
ual intensities do not vary with the intensity of the sources, 
i.e. the fraction of the signal of a source that is subtracted by 
the algorithm is independent of the source brightness. Hence, 
a normalized residual intensity curve obtained by implanting 
artificial point sources of a given brightness can be used to 
calculate the true flux of sources of any brightness and to cor- 
rect the detection limit curve computed from the variance of 
the residual noise. 

The bottom panel of Figure 6 shows that the noise measured 
in the residual image is an adequate measure of the uncer- 
tainty on the intensity of sources at IOct or less. For brighter 
sources (— 25cr), the uncertainty is slightly larger for small 
separations. This is probably due to the larger bias introduced 
by brighter point sources and the more important dependence 
of the amount of partial subtraction on the specific PSF struc- 
ture underlying the point source in regions strongly dominated 
by quasi-static speckle noise. Thus, the noise in the residuals 
may be used as the uncertainty on the flux for most sources 
but it may be necessary to carry out an analysis using artifi- 
cial point sources for brighter sources at small separations. 

4.4. Comparison with previous algorithm 

A comparison of the LOCI algorithm with that used by 
Marois et al. (2006) is presented. Artificial point sources were 
added to the images at several separations in the range 40-500 
pixels (22-275 X/D) by steps of 5 pixels (2.75 X/D). The 
intensities of the artificial sources were adjusted to yield a 
final S/N— 10 with the LOCI algorithm. Both subtraction al- 
gorithms were then run on the images. This was repeated 25 
times with the artificial sources at different azimuthal posi- 
tions. The mean residual intensity and S/N over the 25 az- 
imuthal positions were then computed for each algorithm and 
separation. The results are shown in Figure 7. At all separa- 



6 



Lafreniere et al. 




100 200 300 50 100 150 200 250 

Angular separation {'kID) Angular separation {'kID) 



FIG. 7. — Average normalized residual intensity (top) and S/N ratio (bot- 
tom) as a function of angular separation for the LOCI algorithm (solid line) 
and the algorithm of Marois et al. (2006) (dashed line). 




FIG. 8. — Residual S/N image (including artificial point sources) using 
the algorithm of Marois et al. (2006) (top) and the LOCI algorithm (bottom). 
Both panels are shown with a (-5, +10) intensity range. Each panel is 6"5 by 
3"25. The images have been convolved by a circular aperture of diameter 
equal to W. The saturated region at the center of the PSF is masked out. 

tions, the LOCI algorithm yields a S/N that is better or equal 
to that obtained with the algorithm of Marois et al. (2006). 
The gain is highest at small separations, where it reaches a 
factor ^3, and steadily decreases for larger separations. The 
decrease is most likely due to the increasing relative impor- 
tance of sky background noise. A comparison of the resid- 
ual image of the two algorithms is shown in Figure 8; the 
lower level of noise of the LOCI algorithm is clearly visible. 
The new algorithm yields a better attenuation because it can 
adapt more easily to temporal and spatial variations of the PSF 
quasi-static speckle pattern by using all the images available 
with proper weights (the coefficients) and optimizing the ref- 
erence image combination in smaller subsections. 

The subtraction algorithms were then applied to the original 
sequence of images, i.e. without artificial sources, to compare 
the quasi-static speckle noise attenuation they provide and the 
detection limits they achieve. The quasi-static speckle noise 
attenuation is shown in Figure 9; a single subtraction using the 
LOCI algorithm provides an attenuation of ~10-12 at separa- 
tions of 1-3 arcsec. The formulation of a simple and universal 
criterion for speckle-limited point source detection is usually 



FIG. 9. — Noise attenuation resulting from a single reference image sub- 
traction (bottom) and total noise attenuation (top). The noise attenuation is 
defined as the ratio of the noise in the target image over that in the residual im- 
age; the noise is computed as the standard deviation of the pixel values inside 
an annulus of width ~1 PSF FWHM. The dashed and solid lines are respec- 
tively for the algorithm of Marois et al. (2006) and the LOCI algorithm. The 
attenuations have been corrected for the partial subtraction of point sources. 
Before computation of the initial noise level, a 7 X 7 PSF FWHM median 
filter was subtracted from the images to remove the low spatial frequency 
structures that do not prevent point source detection. 
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FIG. 10. — Statistical distributions of the pixel values of one original S/N 
image after subtraction of a radial profile (dotted line) and of the final S/N 
residual image (solid line) obtained with the LOCI algorithm. From left to 
right, the three panels are for angular separations of 25, 50 and 150 X/D 
respectively. Both images have been convolved by a circular aperture of di- 
ameter equal to W and annuli of area equal to 5000 tt(W/2) 2 were used to 
obtain the distributions at each separation. The continuous solid curve shows 
a Gaussian distribution of unit standard deviation. 

complicated because the distribution of speckle noise is non 
Gaussian (Schneider & Silverstone 2003; Aime & Soummer 
2004; Marois 2004; Fitzgerald & Graham 2006); it possesses 
an important tail at the higher end. However, ADI leads to 
residuals whose distribution closely resembles a Gaussian; 
this is studied in more detail elsewhere (C. Marois et al., in 
preparation). This was indeed verified for the data presented 
here, see Figure 10; a few events above a Gaussian distri- 
bution are seen only at the smallest angular separations. A 
5a threshold is thus adequate for estimating detection limits. 
The final 5a detection limits in difference of magnitudes reach 
13.9, 16.1 and 16.9 at angular separations of 1, 2 and 3 arc- 
sec respectively, see Figure 1 1 . The speckle noise attenuation 
and the detection limits have been properly corrected for the 
partial loss of signal of point sources as measured from the 
residual signal of artificial sources. 

Comparison of the two algorithms were made using a few 
different observation sequences and similar results were ob- 
tained every time. 

5. CONCLUSION 

An algorithm to construct an optimized reference PSF im- 
age used to subtract the speckle noise and improve the sensi- 
tivity to faint companion detection was developed and tested. 
For a given target image limited by speckle noise, the algo- 



PSF subtraction algorithm for high-contrast imaging 
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FIG. 11. — Point source detection limit. The dashed and solid lines are 
respectively for the algorithm of Marois et al. (2006) and the LOCI algo- 
rithm. The detection limits have been corrected for the partial subtraction 
of point sources, for the anisoplanatism observed with ALTAIR and for the 
slight smearing of point sources during an exposure due to FOV rotation. 

rithm linearly combines many reference PSF images such that 
the subtraction of this combination from the target image min- 
imizes the speckle noise. Optimization of the coefficients of 
the linear combination is done for multiple subsections of the 
image independently and the procedure ensures that the min- 
imum residual noise is reached within each subsection. The 
application of the algorithm to ADI yielded a factor of up to 3 
improvement at small separations over the algorithm used in 
Marois et al. (2006). 

The algorithm presented in §2 is general and can be used 
with most high contrast imaging observations aimed at find- 
ing point sources. In particular, it can be used with a sequence 
of images of the target itself obtained at different FOV orien- 
tations (ADI, roll subtraction for HST (Schneider & Silver- 
stone 2003), ground-based observations with discrete instru- 
ment rotations, etc.), with images of the same target at dif- 
ferent wavelengths (simultaneous spectral differential imag- 
ing (SSDI, Racine et al. 1999; Marois et al. 2000) or non- 
simultaneous spectral differential imaging (NSDI) with, for 
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example, a tunable filter) or with images of reference stars 
acquired with the same instrument in a similar configuration. 
The latter could be particularly interesting for HST for which 
the PSF is more stable than at any ground-based telescope 
and for which suitable observations of reference stars may be 
readily retrieved from the archive. This should also be the 
case for the James Webb Space Telescope (JWST), whose 
temperature is expected to be much more stable as a result 
of its more stable environment. Future ground-based instru- 
mentation designed specifically for finding exoplanets will 
have a small FOV, rendering SSDI inefficient to detect plan- 
ets whose spectrum has no steep feature and ADI inefficient 
because of the very long time baseline required for sufficient 
rotation. For such cases, discrete instrument rotations may be 
critical and the algorithm developed here could be used di- 
rectly. The Fine Guidance Sensor onboard JWST (Rowlands 
et al. 2004a), which will include a tunable filter imager (Row- 
lands et al. 2004b) and coronagraph (Doyon et al. 2004), is a 
very interesting prospect for NSDI. Again, the algorithm de- 
veloped here could be applied directly to this case. 
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